GEOG-533 Lab 9

Question 1

Use R package UScensus2010county to complete the following tasks: (20 pt.)

Question 1(a)

Plot a map of New York counties using the plot function.

if (!require(UScensus2010)) {install.packages("UScensus2010")}
library(UScensus2010)
if (!require(UScensus2010county)) {install.county("osx")}
if (!require(UScensus2010tract)) {install.tract("osx")}
library(UScensus2010county)
library(UScensus2010tract)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
df 

plot(shp)

Question 1(b)

Plot a map of New York counties using the qtm function.

library(UScensus2010county)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
df

plot(shp)

if (!require(tmap)) {install.packages("tmap")}
library(tmap)
qtm(shp)

Question 1(c)

How many counties in New York State?

shp <- new_york.county10
df <- shp@data
nrow(df)

[1] 62

Question 1(d)

What’s the 3-digit fips code of Broome County?

shp <- new_york.county10
df <- shp@data
sel <- df$fips[df$NAME10 == "Broome"]
sel

[1] "36007"

sprintf('%03d', as.numeric(sel) %% 1000)

[1] "007"

Question 1(e)

Compute descriptive statistics of the population column (P0010001), including total, minimum, maximum, mean, median, and skewness.

library(UScensus2010county)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
pop <- df$P0010001
sum(pop)

[1] 19378102

min(pop)

[1] 4836

max(pop)

[1] 2504700

mean(pop)

[1] 312550

median(pop)

[1] 91301

library(e1071)
skewness(pop)

[1] 2.517639

Question 1(f)

Create a histogram and a boxplot of the population.

hist(pop)

boxplot(pop)

Question 2

Use R package UScensus2010tract to complete the following tasks: (20 pt.)

if (!require(UScensus2010tract)) {install.tract("osx")}
library(UScensus2010tract)

Question 2(a)

Plot a map of New York census tracts using the plot function.

data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
df 

plot(shp)

Question 2(b)

Compute the total population based on census tracts.

data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
pop <- sum(df$P0010001)
pop

[1] 19378102

Question 2(c)

Select all census tracts in Broome County and plot the map.

data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
sel <- df$county == "007"
plot(shp[sel,])

Question 2(d)

What’s the total population of Broome County?

Broomepop <- df$P0010001[df$county == "007"]
sum(Broomepop)

[1] 200600

Question 2(e)

Create a histogram and a boxplot of population based on census tracts of Broome County.

hist(Broomepop)

boxplot(Broomepop)

Question 2(f)

Select the first five columns of the shapefile of Broome County census tract; add a new population ratio column (= census tract population / county population); save the new shapefile to the hard drive.

dfnew <- df[,1:5]
shp@data <- dfnew
shp@data

shp$ratio <- shp$P0010001/sum(Broomepop)
writeOGR(shp,dsn = "output",layer = "newshp",driver = "ESRI Shapefile",overwrite_layer = TRUE)

NAs introduced by coercion

newshp <- readOGR(dsn = "output", layer = "newshp")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output", layer: "newshp"
with 4919 features
It has 6 fields

newshp@data

Question 3

Use R packages raster and ncdf4 to complete the following tasks: (20 pt.)

if (!require(ncdf4)) {install.packages("ncdf4")}
library(ncdf4)
library(raster)

Question 3(a)

Load the multi-band raster dataset “NDVI.nc” into R.

ndvi <- brick("NDVI.nc")
ndvi

class       : RasterBrick 
dimensions  : 360, 720, 259200, 36  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\NDVI.nc 
names       : X2000.01.01, X2000.02.01, X2000.03.01, X2000.04.01, X2000.05.01, X2000.06.01, X2000.07.01, X2000.08.01, X2000.09.01, X2000.10.01, X2000.11.01, X2000.12.01, X2001.01.01, X2001.02.01, X2001.03.01, ... 
Date        : 2000-01-01, 2002-12-01 (min, max)
varname     : NDVI 

Question 3(b)

Get the basic information about the dataset, including the number of rows, columns, cells, and bands; spatial resolution, extent, bounding box, and projection.

nrow(ndvi)

[1] 360

ncol(ndvi)

[1] 720

ncell(ndvi)

[1] 259200

extent(ndvi)

class       : Extent 
xmin        : -180 
xmax        : 180 
ymin        : -90 
ymax        : 90 

bbox(ndvi)

    min max
s1 -180 180
s2  -90  90

res(ndvi)

[1] 0.5 0.5

projection(ndvi)

[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

Question 3(c)

Aggregate all bands to generate a mean NDVI raster and a maximum NDVI raster; save the two new raster datasets to the hard drive.

ndvimean <- mean(ndvi)
ndvimean

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.8537528  (min, max)

ndvimax <- max(ndvi)
ndvimax

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.9681  (min, max)

writeRaster(ndvimean, "output/ndvimean.nc", overwrite = TRUE) 

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output\ndvimean.nc 
names       : layer 
zvar        : layer 

writeRaster(ndvimax, "output/ndvimax.nc", overwrite = TRUE) 

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output\ndvimax.nc 
names       : layer 
zvar        : layer 

Question 3(d)

Plot the maps of monthly NDVI of the year 2001

ndvi2001 <- ndvi[[13:24]] 
plot(ndvi2001)

Question 3(e)

Create histograms of monthly NDVI of the year 2001.

hist(ndvi2001)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

Question 3(f)

Plot the NDVI map of July 2000; click any location with data on the map and retrieve the NDVI time series for all years; plot the NDVI time series of the selected location.

plot(ndvi, 7)
values <- click(ndvi, n=1, xy=FALSE)

plot(as.vector(values), type = "b")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfError in plot.window(...) : need finite 'xlim' values